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ABSTRACT 

We present simulations of the formation and evolution of galaxy clusters in the Cold 
Dark Matter cosmogony. Clusters with a wide range of mass were selected from pre- 
vious N-body models, and were resimulated at higher resolution using a combined 
N-body/Smooth Particle Hydrodynamics code. The effects of radiative cooling on the 
gas are neglected. While many present-day clusters are predicted to be undergoing 
mergers, the density profiles of those that are approximately in equilibrium are all 
very similar, both for the gas and for the dark matter. These profiles show no sign of 
a uniform density core and steepen gradually from the centre outwards. The standard 
/3-model is a reasonable fit over most of the radius range observable in real clusters. 
However, the value obtained for the slope parameter j3j increases with the outermost 
radius of the fit. Temperature profiles of different simulated clusters are also similar. 
Typically the temperature is almost uniform in the regions which emit most of the 
X-ray fiux but drops at larger radii. The gas temperature and dark matter velocity 
dispersion in equilibrium clusters give values of (3^ = firUpa'jjj^/kT which are con- 
sistent with unity provided an X-ray emission-weighted temperature is used. Larger 
values of i3t are found in merging objects where there is a transient boost in the 
velocity dispersion of the system. Thus /3t > 1 may be an observational indicator 
of merging in real clusters. The similar structure of clusters of differing mass results 
in scaling relations between the X-ray and dynamical properties of clusters identified 
at any given redshift. These scalings are inconsistent with the observed slope of the 
luminosity-temperature relation or the observed sense of evolution of the cluster lumi- 
nosity function. This suggests that the central properties of the intracluster medium 
are determined by non-gravitational processes such as radiative cooling or substantial 
pre-heating at high redshift. 



Oh 

Ql INTRODUCTION 

K/y The dynamical timescales characteristic of galaxy clusters 
are a large fraction of the age of the Universe. As a re- 
sult, cluster formation is an ongoing process which can be 
observed directly in nearby systems, and many properties 
of the initial conditions from which clusters formed are still 
visible in their present structure and dynamical state. De- 
tailed theoretical study of cluster evolution requires mod- 
elling the dynamically dominant dark matter, the hot X-ray 
emitting intracluster gas and the population of cluster galax- 
ies. Although the behaviour of these three components is 
strongly coupled, theoretical studies of galaxy clusters have 
progressed by stages, focussing in the first instance on the 
evolution of the dark matter and on the dynamics of an un- 
evolving galaxy population. Only in the past few years has 
detailed modelling of the dynamics of the intracluster gas 
become possible. 

The simplest model which might capture the essential 
gas dynamics of cluster formation is one in which a nonra- 
diative ideal gas moves in the evolving gravitational field of 
the dissipationless dark matter. This approximation may 



miss some essential physics; although the radiative cooling 
time exceeds the Hubble time for most of the gas in observed 
clusters, much of their X-ray emission often comes from a 
central region where the opposite inequality holds (Fabian 
et al. 1991). In addition, heating and cooling within the 
smaller but denser clumps which were the precursors of ob- 
served clusters could have substantially altered the thermo- 
dynamic state of the gas. Nevertheless, it is clear that a de- 
tailed understanding of the simple problem is a prerequisite 
for studying more realistic (but inevitably more complex) 
situations. This problem also provides a suitable testbed for 
assessing the reliability of the simulation techniques which 
are increasingly being used to investigate cluster evolution. 

A crude but useful analytic model for cluster formation 
assumes spherical symmetry and follows accretion onto an 
initially overdense perturbation. Each shell expands to a 
maximum radius, turns around and falls back onto the per- 
turbation. Shells of dark matter pass through the centre 
and then oscillate within the nonlinear body of the cluster, 
while infalling gas shells encounter an outwardly moving ac- 
cretion shock which thermalizes their kinetic energy. Inside 
this shock the gas is approximately in hydrostatic equilib- 
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rium. Bertschinger (1985) gives a similarity solution for the 
density and temperature structure of the two components 
in the particular version of this infall problem which was 
originally set out by Gunn and Gott (1972). 

If the structure of all clusters is similar, as is expected 
in the above model, then simple scaling relations can be 
derived for the characteristic properties of clusters which 
formed from a scale-free spectrum of initial density fluctu- 
ations (e.g. White 1982). Each cluster present at redshift 
z can be characterised by the mass M contained within a 
sphere of some given overdensity. The mean temperature, 
T, of the cluster and its X-ray luminosity, Lx MpT^^^ 
(p oc (1 -\- z)^ is a characteristic density), then scale as 

T{M,z) X M'-'^^ {1 + z) (1) 
and 

Lx{M,z)ocM'^^^{l + zf^^. (2) 

Kaiser (1986) used these relations to make detailed predic- 
tions for the evolution of the cluster population. 

More recent data have shown that these scaling rela- 
tions fail in a number of important respects. For example 
they imply Lx which is flatter than the observed re- 

lation, Lx y2.7^3.5 (^^(jgg Stewart 1991a, b, Henry 
and Arnaud 1991). Furthermore, at fixed temperature or 
mass the X-ray luminosity is predicted to be larger at earlier 
times, giving rise, for realistic cluster mass functions, to an 
increase in the comoving number density of bright clusters 
with redshift (Kaiser 1991, Blanchard et al. 1992, Bartlett 
and Silk 1993). Current data imply just the opposite trend 
(Edge et al. 1990, Gioia et al. 1990, Henry et al. 1992, Cas- 
tander et al. 1994, Bower et al. 1994). As Kaiser (1991) and 
Evrard and Henry (1991) pointed out, better results are ob- 
tained if the intracluster gas is assumed to have been heated 
to high entropy prior to cluster collapse (perhaps by energy 
injected during galaxy formation) and to have evolved adi- 
abatically thereafter. If the central regions of all clusters 
have the same entropy, then eqn (2) is replaced by 

ix (Af,2) oc M"/^(1 -F^)"/'', (3) 

(Evrard and Henry 1991), so that Lx T^^l^ , independent 
of redshift. 

In recent years several studies have used N-body /hydro- 
dynamic codes (both Lagrangian and Eulerian) to simulate 
the evolution of the intracluster medium in the nonradia- 
tive limit (Evrard 1990a, b; Thomas and Couchman 1992, 
Kang et al. 1994, Bryan et al. 1994). This work relaxed 
the assumption of spherical symmetry and assumed instead 
that clusters form from a Gaussian random field of initial 
density fluctuations. In some cases, the final configurations 
look remarkably like real clusters: the temperature distribu- 
tion is approximately isothermal throughout the body of the 
cluster and the gas density is well fit by the "/9-model" of- 
ten used to characterize observed clusters. However, not all 
simulations give the same results. For example, Thomas and 
Couchman's cluster has a steeper temperature gradient and 
a significantly different density profile than those in Evrard's 
simulations. In addition, there is no consensus on the basic 
issue of whether the simulations obey the expected scaling 
relations. For example, the total X-ray luminosity of the 
rich cluster simulated by Evrard (1990b) seems to increase 



with time, while the cluster population in the simulations of 
Kang et al. (1994) and Bryan et al. (1994) deviates dramat- 
ically from the scaling behaviour predicted in eqns (1) and 
(2). 

In this paper we present new N-body /hydrodynamic 
simulations designed specifically to explore the simple non- 
radiative model, including the case where the gas is initially 
preheated. Our aims are twofold. We wish to test the appli- 
cability of the scaling laws in models with realistic cosmolog- 
ical initial conditions and realistic hydrodynamics. We also 
wish to use the relatively simple behaviour expected theoret- 
ically in this case to assess the reliability and the limitations 
of our numerical techniques. Our initial conditions assume 
a universe dominated by cold dark matter, the paradigm of 
hierarchical clustering theories. We focus primarily on the 
evolution of the distributions of gas and dark matter within 
our clusters. We also examine in some detail the well-known 
"/9-discrepancy" in clusters, and we assess the reliability of 
cluster mass estimates based on X-ray data. Our simulations 
are similar to but have better spatial resolution than those 
presented by Evrard (1990a, b) and Thomas and Couchman 
(1992). 

The rest of this paper is organized as follows. In Sec- 
tion 2 we describe our numerical methods and the param- 
eters of our models. In Section 3 we discuss the evolution 
and the internal structure of our clusters, we compare them 
with Bertschinger's (1985) spherical accretion model, and 
we check how well they obey a similarity scaling. In Sec- 
tion 4 we compare our models with real clusters and we 
discuss the /9-discrepancy. In Section 5 we describe simula- 
tions with pre-heated gas. Finally, we summarize our results 
in Section 6. 



2 THE SIMULATIONS 
2.1 The Code 

Our code is designed to follow the evolution of a mixture of 
coUisional and coUisionless fluids in three dimensions. It is 
based on the Smooth Particle Hydrodynamics (SPH) tech- 
nique and uses a nearest-neighbour binary tree to compute 
gravitational forces. It is fully Lagrangian, free from sym- 
metry restrictions, and has a large dynamic range. It is thus 
well suited to the study of irregular systems with structure 
on a wide range of scales. 

SPH is a Monte Carlo technique which represents a gas 
by "particles" which are used as interpolation centres for the 
computation of thermodynamic quantities and their gradi- 
ents. Each particle is assigned a mass, a density, and a 
temperature; pressure is then obtained from an equation of 
state and gradients are estimated by convolving the parti- 
cle distribution with appropriate kernel functions. In this 
paper we assume an ideal gas, p = (7 — 1)/5m, where p, u, 
and p denote pressure, specific internal energy, and density 
respectively; we take 7 = 5/3. Artificial viscosity terms in 
the equations of motion allow a crude but energetically con- 
sistent representation of shocks. In our code gravitational 
forces between particles are softened using a spline kernel 
of the same form as the interpolation kernel, although the 
corresponding softening length is set independently and is 
kept constant in time. CoUisionless (dark matter) particles 
move under the influence of gravity alone. Excellent reviews 
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of SPH can be found in Benz (1990) and Monaghan (1992), 
a description of the binary tree we use is in Benz et al. 
(1990), our particular implementation of SPH and several 
tests of our code are described by Navarro and White (1993), 
and technical descriptions of SPH as applied to cosmologi- 
cal problems are given by Evrard (1988) and Hernquist and 
Katz (1989). 

2.2 Initial Conditions 

The clusters we have studied were selected from the cos- 
mological N-body simulations carried out by Frenk et al. 
(1990). These simulations modelled cubic regions, 360 Mpc 
on a side, within an Q = 1, A = 0.5 CDM universe. (All 
physical quantities quoted here assume h = 0.5.) Six dark 
matter clumps with one-dimensional velocity dispersions 
ranging from 400 to 1300 km/s were chosen at an output 
time when the rms linear overdensity in a sphere of radius 
16 Mpc took the value, erg = 0.63. The particles in each 
of these clumps were traced back to the initial conditions 
(1 + Zi = 4.74) and a cubic region containing all of them 
was defined. This region was then filled with 10648 new 
particles on a cubic grid which were perturbed using the 
original waves of the Frenk et al. initial conditions, plus 
additional waves chosen to represent the CDM power spec- 
trum between the original resolution limit and the Nyquist 
frequency of the new particle grid. To ensure that each 
cluster in our new simulations is represented by approxi- 
mately the same number of particles, we chose the size of 
this "high-resolution box", L, as indicated in Table 1. The 
distribution of mass outside the high-resolution region was 
modelled by ~ 6000 particles of radially increasing mass us- 
ing a technique similar to that described by Katz and White 
(1993). Overall, this procedure extends the dynamic range 
in mass of the original simulation by a factor of ~ 20 for the 
most massive cluster and ~ 660 for the least massive cluster. 
Evolving our initial conditions with an N-body code we are 
able to recover the clumps in the original low resolution sim- 
ulations, but new small-scale substructure, arising from the 
additional high frequency waves, is sometimes noticeable. 

Within our high resolution region we placed a gas par- 
ticle on top of each of the dark matter particles and gave 
it a very low temperature (~ 15 A ). This ensures that gas 
dynamical effects do not become important until the first 
nonlinear structures form. The mass of a gas particle, nigas, 
is determined by the baryon fraction Qj which we took to be 
0.1. (Note that since cooling is not included, our results can, 
to a good approximation, be directly scaled to other values 
of fij.) The gravitational softening, hg, was set equal to 
100 kpc in all cases. Numerical parameters for all our sim- 
ulations are listed in Table 1. The symbols will be defined 
in the text when they are first mentioned. 

In addition to our six standard runs, labelled Cil-6, 
we ran three further simulations in which the gas was "pre- 
heated" at 2 = 3. This was accomplished by raising the 
temperature of each gas particle so as to give a uniform 
specific entropy, s = ln{T/ p"^^^), equal to the specific en- 
tropy at the centre of CLl, our most massive cluster, at 
z = 0. Our three "preheated" models are labelled EPl to 
EPS in Table 1 and have the same initial conditions as the 
standard models CLl, CL3, and CL5 respectively. Simu- 
lations CLl X 3 and Cil/2 use the same initial conditions 



as that of run CLl, but the total number of particles has 
been increased by a factor of ~ 3 and decreased by half, 
respectively. 

2.3 Estimates of X-ray luminosity and code 
units 

The bolometric X-ray luminosity of a cluster is given by the 
volume integral of the emissivity. 



Lx = / L{T)dV, 
'v 



(4) 



where L(T) = ra^A(T). H ere, the gas number density, n = 
p/finip, where nip is the proton mass, fi = 0.6 for a fully 
ionized primordial plasma (Spitzer 1978), and A(T) is the 
cooling function which, at the temperatures relevant to rich 
clusters, can be well approximated by 

1/2 



A(T) = 1.2 X 10' 



-24 



T 



3 -1 
erg cm s , 



(5) 



erg s 



-1 



(6) 



IkeV^ 

For the simulations we use the estimator, 
Lx = 1.2 X 10'^"^ {fimp)~^mg ^ ..-^^^^ 

8 = 1 

where nig is the mass of a gas particle, and pi and Ti are 
the density and temperature (in keV) at the position of the 
i-th gas particle. 

The computational units of the code are given by: 



G 



[mass 
[distance 

[time 

[velocity 
[density 
[energy /mass 

[energy /time 



1 

10^° m„ Mq 
1 du kpc 



: 4.71 X 10° 



207.4 



du 



1/2 



1/2 



yrs 



km/s 



6.77 X 10 I J g/cm 



:4.3xl0"f^j erg/g 
X 5/2 

5.75 X 10 I — — I erg/s. 
du j 



Physical quantities from the simulations can be consistently 
rescaled to other situations by changing the values of rau and 
du- The numbers quoted throughout this paper correspond 
to rau = du = 1- 



RESULTS 



3.1 Evolution 

Figure 1 shows the evolution of the distributions of dark 
matter and of gas in CLl, our most massive cluster. Over- 
all, the two distributions look very similar although, as dis- 
cussed below, there are significant differences of detail. The 
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Table 1. Cluster properties (2: = 0). 



Model 


L 


M200 




^200 


^ DM 


gas 


Rons / Rr> A/f 


Tx 


Lx 




[Mpc] 


L \y A 


[lol^Mo] 

L \:j i 


[Mpc] 


[km/s] 






[K] 


[erg/ s] 


CLl 


45.0 


2.29e5 


5.94 


3.43 


1.17e3 


0.093 


1.17 


1.10e8 


3.98e45 


CL2 


38.2 


2.19e5 


3.58 


3.36 


1.15e3 


0.087 


1.11 


8.13e7 


2.34e45 


C'L3 


31.0 


1.05e5 


1.94 


2.63 


8.91e2 


0.086 


1.11 


4.57e7 


1.02e45 


n T A 




. ioe^y: 


11/1 


Z .Uo 




U.Uo / 


1 91 




. / oe^yz^y: 


CL5 


20.9 


1.15e4 


0.59 


1.26 


4.37e2 


0.091 


1.09 


1.26e7 


6.46e43 


CL6 


14.0 


9.77e3 


0.18 


1.19 


3.72e2 


0.089 


1.22 


1.00e7 


4.07e43 


CLl X 3 


45.0 


2.34e5 


1.93 


3.44 


1.17e3 


0.097 


1.21 


l.lle8 


4.36e45 


CLl/2 


45.0 


2.43e5 


11.9 


3.48 


1.18e3 


0.094 


1.33 


1.15e8 


2.88e45 


EPl 


45.0 


2.29e5 


5.94 


3.43 


1.17e3 


0.094 


1.19 


1.17e8 


3.09e45 


EP2 


31.0 


1.05e5 


1.94 


2.63 


8.91e2 


0.085 


1.17 


6.02e7 


6.16e44 


EP3 


20.9 


1.12e4 


0.59 


1.26 


4.27e2 


0.071 


1.67 


1.86e7 


1.51e43 



cluster grows primarily through mergers, as predicted by the 
analytic model originally due to Press and Schechter (1974), 
and extended by Bond et al. (1991) and Bower (1991). At 
z = 0.84 and 0.40, the fraction of the present mass of the 
cluster that is in clumps (identified using a friends-of-friends 
algorithm) of more than ~ 100 particles is 0.48 and 0.75, re- 
spectively, compared with 0.37 and 0.68 predicted by the 
analytic Press-Schechter theory. Increasing the number of 
particles of the simulation by a factor of three (run CLl x 3) 
brings these numbers down to 0.45 and 0.73, improving the 
agreement with the analytic model. At each stage gas col- 
lapses into clumps with the dark matter and is shock-heated 
to the relevant virial temperature. This hierarchical growth 
continues to the present and results in significant irregularity 
in our final clusters. Figure 2 shows the final configuration 
of the gas in clusters C L2-C L6. In several cases, substruc- 
ture is still clearly discernible and some of the clusters have 
massive companions nearby. 

The rate at which mass is added to the clusters is il- 
lustrated in Figure 3. We plot the evolution of the dark 
mass within a sphere of current mean overdensity of 200. At 
each time we centre this sphere on the most massive clump 
present. The curves are normalized to the mass M'joo con- 
tained at the final time, and different symbols correspond 
to different simulations. Mergers stand out in this plot as 
sudden jumps in mass. As expected in an Q = 1 universe, 
the formation of rich clusters is a very recent phenomenon. 
Half the clusters in Figure 3 have doubled their mass since 
z = 0.2, and most have undergone substantial mergers since 
that time. This explains the substantial irregularity seen in 
Fig. 2. As pointed out by Forman and Jones (1990), a large 
fraction of real X-ray clusters show significant substructure. 
Several authors (eg Richstone et al. 1992, Lacey and Cole 
1993, Kauffmann and White 1993, Evrard et al. 1994b) have 
noted that this argues in favour of a high value of Qq- 

Figure 4 shows the evolution of the thermodynamic 
properties and the X-ray luminosity of the intracluster gas. 
The gas which ends up near the median radius of a clus- 
ter typically reaches high density as the first resolved struc- 
tures form; its average density remains essentially constant 
thereafter (dotted lines). The same is almost true for the 
gas which ends up in a cluster core; its mean density in- 
creases only slightly with time after the initial phase of evo- 



lution (solid lines). The mean density of the gas which ends 
up in the outer regions actually declines slowly over most 
of the evolution (dashed lines). The mean temperature of 
the gas which ends up in each of the three regions shows a 
slow but steady increase due to adiabatic compression and 
shock-heating. The shocks caused by accretion and mergers 
tend to be mild because the relative speeds are generally 
not highly supersonic. This is particularly so in the central 
regions where the entropy increases only slowly after the 
collapse of the first subclumps. In the outer parts shocks 
due to recent mergers are more effective in raising the en- 
tropy. The X-ray luminosity is determined primarily by the 
gas density near the core and depends only weakly on tem- 
perature (eqs 4-6). On average it changes little after z = 1, 
but it fluctuates whenever a merger drives oscillations in 
the central density. (These can have amplitudes of up to 
a factor ~ 3.) After a merger the X-ray luminosity is gen- 
erally slightly higher than the combined luminosity of the 
progenitors. 

The results shown in Figure 4 are similar to those of pre- 
vious simulations (see e.g. Figure 1 of Thomas and Couch- 
man 1992 and Figure 4 of Evrard 1990b). Notice, however, 
that they do not seem consistent with the predictions of 
the simple scaling laws we gave in the Introduction. For 
example, these predict the luminosities of clusters to de- 
crease with time. We now investigate in detail the internal 
structure of our clusters and compare their properties and 
evolution with the predictions of the spherical infall model 
and the self-similar scaling relations. 

3.2 Density profiles 

At the present epoch several of our clusters have substantial 
substructure within the virial radius which reflects recent 
mergers (see Figure 2). In the present section we wish to 
study the "equilibrium" state of clusters. Hence for those 
clusters with significant substructure at 2 = we consider 
slightly earlier epochs when the most recently accreted mas- 
sive clump was still outside the virial radius. In all cases we 
can choose such a time at z < 0.2. The results presented in 
this and the following two subsections refer to these "equi- 
librium" configurations. It is important to note that they 
cannot he considered typical of nearby clusters. Rather they 
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Table 2. Density profile fits. 



Model 


Ao 


' c 




r 


Ao 


' c 


% 




[gas] 


[gas] 


[gas] 


[gas] 


[DM] 


[DM] 


[DM] 


CLl 


1.57e3 


0.106 


0.86 


1.30 


2.40e4 


0.082 


0.86 


CL2 


2.77e3 


0.101 


0.77 


1.40 


9.17e4 


0.054 


0.75 


CL3 


4.22e3 


0.101 


0.90 


1 38 


2 91e4 


125 


1 05 


CL4 


1.89e3 


0.074 


0.74 


1.35 


3.20e4 


0.066 


0.78 


CL5 


1.28e3 


0.127 


0.84 


1.15 


2.42e4 


0.163 


0.90 


CL6 


1.31e3 


0.098 


0.80 


1.42 


1.76e4 


0.093 


0.85 


' rc is in units of 


^200- fits are 


within 


'■200- 







are typical of nearby regular clusters, which may be a mi- 
nority of all clusters. 

Density profiles averaged over spherical shells are shown 
for the gas and the dark matter in Figures 5a and 5b respec- 
tively. The center of each system was defined by an iterative 
procedure that computes the center of mass inside a spheri- 
cal region and then shrinks that spherical region until about 
100 particles are left. This picks up the center of the most 
massive clump in the volume of interest. The two compo- 
nents have similar profiles which are poorly described by 
a single power-law. Their logarithmic slope steepens from 
~ —1 near the centre to ~ —3 at large radius. There is no 
evidence for a constant density core outside the region where 
the gravitational softening is important. (At twice the soft- 
ening radius the density in all clusters has dropped by at 
least a factor of two from its value at hg.) Scaled density 
profiles for the six clusters are shown in Figure 6a. Here we 
plot overdensity (density divided by the mean density of the 
universe at the time when the cluster is identified) as a func- 
tion of the dimensionless radial variable, r/vjoo, where r2oo 
is the radius of the sphere within which the mean interior 
overdensity is 200. The scaled profiles are all remarkably 
similar outside the region affected by gravitational soften- 
ing. Bumps in the profiles beyond r2oo are due to compan- 
ion systems which are falling onto the cluster. Despite the 
lack of a constant density core, the scaled density profiles of 
both gas and dark matter are fitted moderately well by the 
"/?-model". 



p/7i = Ao{l + (r/rcf 



-3/3//2 



(7) 



over all but the innermost regions. This fit is shown as the 
thick solid line in Figure 6b which has = 0.8, Ao = 
1.4e3 and Tc/t'joo = 0.1 for the gas and = 0.8, Ao = 
2.2e4 and Tc/t'joo = 0.09 for the dark matter. The best- 
fit parameter values for the individual clusters are listed in 
Table 2. (Note that while the slope of the outer profiles is 
quite well determined, the values of Tc and Ao are affected 
by poor resolution near the centre and have highly correlated 
uncertainties.) 

While the /9-model has traditionally been used to de- 
scribe cluster density profiles, it has little theoretical jus- 
tification and is not a very good fit to the inner regions, 
particularly for the dark matter. A better fit to the dark 
matter profile for r < r2oo is given by the simple fitting 
function 



P/P 



1500 4oo 
r (5 r + T'joo)^ 



(8) 



This function behaves as at small radius in agreement 
with our best resolved models and with the much higher 
resolution simulations of Dubinski and Carlberg (1991). 

The thick dotted lines in Figure 6b show Bertschinger's 
(1985) self-similar infall solution. This describes spherical 
accretion of coUisionless and coUisional fluids onto a point 
mass in an otherwise unperturbed Einstein-de Sitter uni- 
verse. It predicts p oc r~^^^ at small radius. The simulated 
profiles are somewhat shallower than the similarity solution, 
particularly for the dark matter. This is not completely 
surprising, for profiles shallower than p oc r~^^^ have al- 
ready been found in several analytic and numerical works 
(Fillmore and Goldreich 1984, Quinn et al. 1986, Hoffman 
1988, and Efstathiou et al. 1988). The discrepancy is most 
marked near the centre where the simulated profiles fall well 
below Bertschinger's solution. This appears to be a genuine 
disagreement since it persists to radii several times the soft- 
ening length in our best resolved clusters. However, still 
higher resolution simulations are needed for a definitive de- 
termination of central cluster profiles. 

The scaled density profiles of the largest progenitors of 
the C Ll-C L6 clusters at 2 = 1 are shown in Figure 7. Be- 
yond the region dominated by the gravitational softening, 
these profiles are very similar to those at 2 = 0, which are 
illustrated by the thick solid lines. (Note that the slopes 
are the same in these regions; the small offset is due to our 
definition of r200) which forces all the curves to encompass 
the same mass within this radius.) Near the center, the sim- 
ilarity breaks down, due to the smaller number of particles 
per cluster (up to a factor of ten), and to the fact that the 
gravitational softening is now a much larger fraction of r2oo- 
As a result, the central overdensities are on average about 
three times lower than at 2 = 0. This can lead to systematic 
underestimation of physical quantities that depend crucially 
on the central properties of the clusters, such as the X-ray 
luminosity (see §3.6 below). 

Although the gas and dark matter profiles appear very 
similar in Figure 6, the gas is actually slightly more extended 
than the dark matter. Inside r2oo the gas fraction ilgas is 
always smaller than the global value of 0.1 (see Table 1). 
This difference is also apparent in the ratio of the half mass 
radii of the gas and dark matter, Rgas / Rdm t which, as 
noted in Table 1, is typically greater than unity. (Note that 
these radii are determined relative to the total mass within 
''200-) This property reinforces the conclusion of White et al. 
(1993) that the large gas fractions observed in real clusters 
cannot be reconciled with the much smaller global fraction 
expected from Big Bang nucleosynthesis considerations in 
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an Q = 1 universe. 

The segregation between dark matter and gas results 
from the different behaviour of the two components during 
a merger. Shocks prevent colliding gas distributions from 
interpenetrating. This produces a transient offset between 
gas and dark matter and a net transfer of energy and angu- 
lar momentum to the gas (Navarro and White f993). The 
effect can be illustrated by following the mean properties of 
the particles that end up in the core of a cluster, defined to 
be the region containing the inner fO% of its mass. The tra- 
jectory of each cluster core from 2 = f to 2 = is shown in 
the density-temperature or density-velocity dispersion plane 
in Figure 8. The upper dotted line shows the track corre- 
sponding to isentropic gas, pgas T^/^^ while the lower 
dotted line corresponds to constant phase-space density for 
the dark matter, Pjjm f^. Shocks increase the entropy of 
the core gas and cause temperatures to rise with relatively 
little change in density. In contrast, the dark matter evolves 
at roughly constant phase-space density. As a result, the 
gas density in the cluster core is only 0.07 times the dark 
matter density, whereas a ratio ofQg/(f — Q^) = O.ff would 
be expected if the two distributions were identical. 

3.3 Temperature profiles 

Figure 9 shows scaled gas temperature profiles for our six 
clusters. These were obtained by normalizing to the charac- 
teristic temperature, T200 = P-mpGM20ol {2kr20o), where k 
is Boltzmann's constant. The profiles agree to within about 
a factor 2, a similar level to that seen earlier for the density 
profiles. (Note the very different scales in the ordinates of 
Figures 6 and 9.) All the clusters have a near-isothermal in- 
ner region that extends, in the best resolved case, for about a 
decade in radius. Near r ~ 0.4r200) the temperature begins 
to drop and by r2oo it has fallen to about half the central 
value. This is as expected for gas in hydrostatic equilib- 
rium given the density profiles of Figure 6. In the region 
0-1 ^ i'/i'200 ^ 0-4, the gas and dark matter density pro- 
files have an effective slope of ~ — 2 (indicated by the thick 
solid line in Figure 6a), whereas near r/vjoo — f , the profiles 
steepen to r~^'^. The equilibrium temperature is therefore 
nearly constant in the inner region and drops approximately 
as T~^l'^ in the outer region. 

The thick dotted line in Figure 9 shows the gas tem- 
perature in Bertschinger's (f 985) similarity solution. In this 
model the temperature increases inwards from the shock and 
approaches an r~^l^ power-law near the centre. In the outer 
cluster this model gives an acceptable match to the simula- 
tions, but in the inner parts it overestimates the temperature 
just as it overestimates the density (see Figure 6). 

The temperature profiles that we find agree with that 
of Evrard (f990b) but disagree with that of Thomas and 
Couchman (f992). (Since the average temperature of the 
cluster in each of these studies is ~ f 0®K, the corresponding 
virial radius is r2oo — 3 Mpc.) In Evrard's calculation the 
near-isothermal region extends almost to r200, as is the case 
in one of our own models, whereas the temperature drops by 
almost an order of magnitude between 0.fr2oo and r2oo in 
Thomas and Couchman's cluster. None of our clusters shows 
such a dramatic decline in temperature. The reason for this 
discrepancy is unclear. All three studies assume the same 
cosmology and attempt to model the same physics using 



similar codes. Both Evrard and Thomas & Couchman select 
clusters as high peaks of the initial linear density field while 
we select clusters from the final nonlinear mass distribution. 
This seems unlikely to be the source of the problem given 
the dramatic difference between the two earlier simulations. 
It is possible that Thomas and Couchman's cluster is simply 
an atypical example. 

3.4 Hydrostatic equilibrium and binding mass 
estimates 

Using the density and temperature profiles derived in the 
preceeding subsections we can quantify how well hydrostatic 
equilibrium describes the intracluster medium of our simu- 
lated clusters. For real clusters the interpretation of X-ray 
data, and in particular the estimation of cluster mass, is 
usually based on the assumption that the gas is in equilib- 
rium. Our simulations allow us to determine the accuracy 
of such estimates and to assess the biases which may arise 
when only incomplete information is available. 

For a spherically symmetric cluster with the density 
profile of eqn (7), the equation of hydrostatic equilibrium 
may be written as 

MTir) = - I) ' ' +- T r r. 9 

In Figure f we compare the mass derived from this equation 
with the actual mass, Mtot, contained within radius r. Out- 
side the region where gravitational softening is important, 
eqn (9) underestimates the true mass by up to ~ 30%. As 
pointed out by Evrard (f990b), this discrepancy reflects the 
existence of bulk motions in the gas whose kinetic pressure 
makes a non-negligible contribution to the total pressure 
support. 

The importance of kinetic pressure in the simulations 
may be quantified by the parameter jSgas = t^'n^pf^gas/kT, 
where agas is the velocity dispersion of the gas particles. 
For hydrostatic equilibrium, Pgas = 0. The values of Pgas 
for our simulated clusters (Table 3) range between ~ O.f-0.3 
and are large enough to account for the discrepancy between 
the true mass and that estimated from eqn (9) and for the 
fact that the measured temperatures in Figure 9 fall slightly 
below the "virial temperature." The magnitude of the resid- 
ual bulk motions depends on the detailed dynamical history 
of the system so a constant "fudge factor" for the mass es- 
timates, as proposed by Evrard (f990b), is a rather crude 
approximation. Systems that have not experienced recent 
mergers (eg the solid curve in Figure fO) satisfy eqn (9) 
very well at all radii. 

Although the ASCA X-ray satellite is now beginning 
to measure gas temperature profiles in clusters, for most 
objects only a single emission-weighted temperature is cur- 
rently available. In this case, the gas is commonly assumed 
to be isothermal and the logarithmic derivative of the tem- 
perature in eqn (9) set equal to zero. The resulting mass 
estimates are shown in the lower panel of Figure f 0. In most 
cases, the true mass is underestimated by a similar factor 
to that obtained when the temperature profile is known, re- 
flecting the fact that the temperature distribution of our 
clusters is close to isothermal in the region where most of 
the X-ray emission is produced. 
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3.5 Scaling Laws 

We now examine how well our simulated clusters obey the 
scaling relations described in the Introduction. The filled 
circles in Figure 11 show correlations between X-ray and 
dynamical properties for clusters C Ll-C L6. (The filled tri- 
angles correspond to our "preheated" gas simulations and 
the crosses to data for real clusters; these will be discussed 
below.) Here, we consider the properties of the clusters at 
z = 0, rather than those of the "equilibrium" configurations 
discussed in the three preceeding subsections. Numerical 
values of the quantities plotted are given in Table 1. 

The solid lines in Figure 11 show the scaling relations 
derived from eqns (1) and (2). The mass-temperature re- 
lation is normalized assuming that the mean overdensity 
of the clusters is 200: T = 5.1(M2oo/lO^®M0)2/^keV. 
The mass-velocity dispersion is normalized assuming that 
the clusters are isotropic isothermal spheres: iJoM = 
(0.5GM2oo/''20o)"'^''^ km/s. The calculation of the X-ray lu- 
minosity requires additional assumptions about the struc- 
ture of the clusters (e.g. the size of the core radius), so we 
fix the normalization of the mass-luminosity relation empir- 
ically by fitting to the results of the simulations. This gives 
Lx = 1-2 X lO'*^(M2oo/lO^®M0)'*/^erg/s. The normaliza- 
tions of the other scaling relations follow from these. Over- 
all, our simulations follow the scaling relations remarkably 
well. There is a small offset in the T-M'joo relation which 
results from the fact that the gas is incompletely thermal- 
ized and thus not in strict hydrostatic equilibrium. This also 
affects the T-udm relation. 

3.6 Sensitivity to numerical parameters 

The resolution of our simulations is determined primarily 
by the total number of particles, N , which effectively sets 
the minimum value of the gravitational softening, hg, re- 
quired to suppress two-body relaxation. Figure 12 shows 
the effect of increasing and decreasing the number of par- 
ticles in the CLl cluster by factors of 3 and 2 respectively. 
The largest simulation (dotted lines, run CLl x 3) has total 
N ~ 70, 000 and the smallest one (dashed lines, run Cil/2) 
has N ~ 11, 000. The solid lines correspond to the original 
CLl simulation. Even though the total number of particles 
varies by a factor 6, the density and temperature profiles are 
practically indistinguishable. 

As shown in Tables 1 and 3, dynamical properties are 
also in good agreement in the three runs (CLl, CLl x 3 
and Cil/2). For example, the residual bulk motions of the 
gas iPgas in Table 3) are fairly insensitive to N and thus 
seem to be unaffected by discreteness effects in the evolu- 
tion of the cluster potential. Even the X-ray luminosity, 
the quantity most sensitive to numerical resolution because 
of its strong dependence on central density, seems to have 
converged. Increasing TV by a factor 3 increases Lx by less 
than 10% although decreasing TV by a factor of 2 causes Lx 
to drop by ~ 30%. Decreasing the number of particles even 
further results in severe underestimates of the total X-ray 
luminosity of the cluster. This is shown in Figure 13, where 
we plot how Lx depends on the total number of gas par- 
ticles inside T'joo (Ngas)- Only when Ngas is larger than 
about 2, 000-3, 000 do the X-ray luminosity estimates seem 
to be reliable. 



Other factors that may influence the central regions 
of a cluster, and therefore might spuriously affect the X- 
ray luminosity of a simulated cluster, include the choices of 
gravitational softening, initial redshift, and SPH smoothing 
lengths. We have carried out a number of experiments to 
test these possibilities. Varying hg from 100 to 200kpc, 1 + Zi 
from 4.74 to 9.48, and requiring 100 neighbours instead of 
the usual 32 in the SPH smoothing algorithm cause, how- 
ever, little change in the X-ray luminosities at 2 = for 
simulations with our standard particle number. Therefore, 
we conclude that the density profiles of these clusters have 
a well-defined "core" that is resolved in our simulations at 
z = 0. (We use the term "core" here rather loosely to sig- 
nify the radius where the total X-ray luminosity starts to 
converge.) Inside this "core" the density is not constant, 
but its logarithmic slope is relatively shallow, so that once 
the numerical resolution is good enough, further improve- 
ments have no appreciable effect on Lx ■ Note that were this 
"core" produced by the gravitational softening we would not 
recover the Lx oc relation predicted by the scaling laws, 
since hg does not scale from cluster to cluster. The gravita- 
tional softening was kept constant in physical units, and so 
the softened region is a different fraction of the virial radius 
in different clusters (cf. Figure 6). 

The effects of limited numerical resolution can have 
dramatic effects on the evolution of the X-ray luminosity. 
At early times clusters are less massive and therefore more 
poorly resolved than at late times. This problem is com- 
pounded by the fact that at early times the gravitational 
softening, which is fixed in physical units, is a larger fraction 
of the virial radius. As shown by the z = 1 density profiles 
in Figure 7, these two effects significantly limit the central 
density that a cluster can attain. At this redshift our clus- 
ters contain typically no more than ~ 400 particles inside 
their virial radii. Figure 13 shows that this must severely 
compromise numerical estimates of their X-ray luminosity of 
a cluster, leading to a systematic underestimation of Lx at 
high redshift and to a breakdown of the scaling behaviour. 
An effect of this magnitude is expected due to lack of res- 
olution at this redshift (see Figure 13). This is illustrated 
in Figure 11 where we plot the Lx-M'joo relation for the 
progenitors of clusters C Ll-C L6 aX z = 1 (open circles). 
The dashed line shows the expected behaviour predicted by 
eqn (2). The redshift evolution of the simulated clusters 
is much weaker than the (1 -|- z)^ 1"^ dependence implied by 
the scaling laws; they are about 4 times less luminous than 
predicted. An effect of this magnitude can be entirely due 
to lack of resolution at this redshift, as shown in Figure 13. 
Although based on the present set of simulations we cannot 
demonstrate that the redshift evolution agrees with the the- 
oretical expectations, the similarity of the "scaled" clusters 
in the regions not affected by numerical limitations leads us 
to suggest that the profiles shown in Figure 6 are generic for 
clusters formed in this model, and that the departure of sim- 
ulated clusters from the predicted evolutionary behaviour is 
likely to be due mainly to resolution effects. 

The dependence of the estimated X-ray luminosity on 
numerical resolution is not surprising but has sometimes 
been overlooked. Kang et al. (1994) and Bryan et al. (1994) 
interpreted a discrepancy between the redshift evolution of 
the X-ray luminosity and the scaling laws as a genuine pre- 
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diction of their models rather than as a reflection of their nu- 
merical limitations. For the same reason, their conclusions 
regarding clusters at the present epoch may also be affected 
by resolution effects. Since their clusters were identified in 
single simulations of cosmological volumes, they are resolved 
to varying degrees. Poorer clusters have smaller numbers of 
particles and subtend fewer resolution lengths. As a result 
their central densities and X-ray luminosities may be sys- 
tematically underestimated by a larger factor than those of 
richer clusters. This would lead to a steeper Lx-T relation 
than the scaling laws predict. At 2 = 0, our simulations are 
not affected in this way because we use approximately the 
same number of particles to simulate clusters of differing 
mass. Our derived Lx-T relation thus agrees remarkably 
well with the theoretical scaling relation. 

4 COMPARISON WITH REAL CLUSTERS 



4.1 Gas density profiles 

The gas density profiles discussed in Section 4.2 are mod- 
erately well fit by the "/9-models" (eqn 7) often used to 
describe observed clusters. Since the core regions of real 
clusters are affected by cooling processes which we have ne- 
glected, we limit our comparison to the outer profile, as mea- 
sured by the parameter [ij. As noted by Evrard (1990a), 
typical values of [ij found in simulations of this kind are 
~ 0.8, rather than ~ 0.6, the value inferred for most ob- 
served clusters from their X-ray surface brightness profiles. 
However, this comparison can be misleading. Since eqn (7) 
is not a perfect description of our models, the "best" value 
of [ij depends on the range of radii used for the fit. We illus- 
trate this in Figure 14 by plotting the value of [ij derived for 
the ensemble of Figure 6b against Tmax, the outermost ra- 
dius used in the fit. While we find P ^ ~ 0.8 for Tmax — T'joo, 
/]f is only ~ 0.6 for Tmax — 0.2r2oo 

The value of Tmax used in analyses of real clusters de- 
pends on the quality of the data. The filled circles in Fig- 
ure 14 show values oi P ^ for those clusters from the compi- 
lation of Jones and Forman (1984) for which we have been 
able to infer the value of Tmax- We assume Tmax to be 
the radius at which uncertainties in the measured surface 
brightness exceed a factor of 10 (see Figures 1 and 2 of 
Jones and Forman; we take r^oo ~ 3(Tx/7.5keV)"'^''^ Mpc 
for their clusters). The point at Tmax /t'joo = 1 corresponds 
to the ROSAT observations of the Coma cluster reported 
by Briel et al. (1992). The trend in the real data is similar 
to that found in the simulations. In most cases, the obser- 
vations become noisy well inside the virial radius and this 
may account for the low [ij values measured. As discussed 
by Evrard (1990a), a further bias could result from incorrect 
background substraction. The [ij values for real clusters lie 
slightly below those for the simulations but overall there is 
reasonable agreement. Higher signal-to-noise observations 
should help to clarify these issues. 

4.2 Temperature profiles 

The temperature profiles of our clusters are characterized 
by a near-isothermal region extending to r ~ 0.5r2oo and 



a gradual decline to about half the central value near r2oo- 
Profiles for many real clusters should soon be available from 
observations with the ASCA satellite, but for the time being 
we compare our results with the Coma cluster for which a 
radial temperature profile has been measured to ~ 1 Mpc 
(Watt et al. 1992), about 1/3 of the virial radius. In this 
region, the cluster is very close to isothermal. Although the 
temperature further out has not been directly measured, 
detailed modelling of X-ray and optical data by Hughes et 
al. (1988a, b) suggests a drop in temperature similar to that 
seen in our models. Hughes et al. derive a polytropic index, 
r ~ 1.6 for the outer parts of Coma. Values of F obtained 
by fitting a power law to the temperature profiles of our 
clusters in the range 0.4 < r/vjoo < 1 (see Figure 9) are 
listed in Table 2. These values are somewhat lower than the 
value quoted by Hughes et al. but, given the considerable 
uncertainty in the analysis of the data, we conclude that, 
as with the density profiles, the temperature profiles of our 
simulated clusters are consistent with those inferred for real 
clusters. 

4.3 Correlations 

X-ray luminosities, gas temperatures and galaxy velocity 
dispersions for various samples of clusters are plotted in the 
correlation diagrams of Figure 11. The crosses represent 
data from Edge and Stewart (1991a, b), Henry and Arnaud 
(1991), David et al. (1993), and Lubin and Bahcall (1993). 
Starred polygons correspond to two Hickson groups (HG92 
and HG62) recently detected with ROSAT (Mulchaey et al. 
1993, Ponman and Bertram 1993). 

The fundamental Lx-T relation for the data obeys 
the well known law, Lx9ST^, with the two Hickson groups 
falling on the extrapolation of the relation for rich clusters. 
This is significantly steeper than the corresponding relation 
for our simulated clusters which closely follow the scaling 
law of eqn (2), Lx oc T^ . Thus, at low temperatures, our 
models give about an order of magnitude higher X-ray lu- 
minosity than is observed. The discrepancy between our 
models and the observed Lx-T relation is reflected in the 
Lx-<y plane although in this case the scatter in the data 
is larger. By contrast, the observed T-a relation does not 
show any systematic deviation from the relation obeyed by 
our models, T oc , appropriate to an isothermal gas in 
hydrostatic equilibrium. However, the scatter in the data 
is larger than in the models, a feature related to the "/?- 
discrepancy" which we discuss in some detail below. 

Several suggestions have been put forward to explain 
why the simple nonradiative models fail to account for the 
observed Lx-T relation. One possibility is that the gas frac- 
tion increases with cluster temperature (David et al. 1990, 
Arnaud et al. 1992). Another is that cooling preferentially 
depresses the luminosity of poorer clusters. A third (which 
we discuss further below) is that the central gas densities 
and core radii vary systematically with temperature, as ex- 
pected if the intracluster medium were heated to relatively 
high entropy at early times. 

4.4 The "/3-discrepancy" 

The "/9-discrepancy" in clusters is a problem which has ex- 
ercised X-ray astronomers for many years. It refers to a dis- 
crepancy which arises when attempting to fit cluster data 
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Table 3. Values of the jSj' parameter. 



Model 


Pt(Tx) 






Pt(Tx) 


/?T (^m) 


f^aas {Tx ) 


Pt(Tx) 




[''200] 


[''200] 


[''200] 


[''1000] 


[''1000] 


[''1000] 


(z = 0) 


CLl 


1.01 


1.18 


0.11 


1.08 


1.10 


0.07 


1.01 


CL2 


1.07 


1.14 


0.15 


1.15 


1.17 


0.11 


1.28 


CL3 


1.18 


1.36 


0.36 


1.31 


1.32 


0.37 


1.37 


CL4 


1.04 


1.18 


0.22 


1.08 


1.18 


0.12 


1.24 


CL5 


1.05 


1.23 


0.08 


1.24 


1.16 


0.06 


1.18 


CL6 


1.09 


1.38 


0.30 


1.13 


1.10 


0.19 


1.09 


CLl X 3 


0.94 


1.17 


0.10 


1.05 


1.11 


0.06 


0.94 


CLl/2 


0.92 


1.12 


0.09 


1.00 


1.06 


0.07 


0.92 



to the /9-model proposed by Cavaliere and Fusco-Femiano 
(1976). This model assumes that clusters are spherical, 
are in hydrostatic equilibrium, and that both the intraclus- 
ter medium and the galaxies may be modelled as isotropic 
isothermal "gases" with density profiles given by eqn (7). 
With these assumptions the hydrostatic equilibrium equa- 
tion can be integrated immediately to give the gravitational 
potential of a cluster. 



$(r) 



2 /imp 



ln(l + r Itc) 



(10) 



Note that the density corresponding to this potential falls 
as r~'^ at large radii and that it cannothe cast in the form 
of equation (7). Since gas and galaxies are in equilibrium 
in the same potential well, equation (10) shows that the 
two components must have the same value of the product 
l^jT. Thus if we retain [ij and T for the properties of the 
gas, and denote the corresponding properties of the galaxy 
distribution by f)j^, and T* = /impa^/k where fr* is the one- 
dimensional velocity dispersion, the slopes of the density 
profiles of the two components are related to the ratio of 
their temperatures by 



kT 



(11) 



The galaxy distribution is often modelled by King's approx- 
imation to the core of an isothermal sphere (eqn 7 with 

= 1). In this case, the slope of the gas density profile, 
Pf, should equal Pj'. As we saw earlier, fits to X-ray sur- 
face brightness profiles typically give ~ 0.6 (Bahcall and 
Sarazin 1977, Mushotzky 1984, Jones and Forman 1984), 
whereas measurements of gas temperatures and galaxy ve- 
locity dispersions give Pj' ~ 1. (Many values of Pj' which 
were found to be significantly greater than one in early stud- 
ies have subsequently been revised downwards as a result 
of improved velocity dispersion determinations. Edge and 
Stewart (1991a,b) obtained = 0.91 ± 0.08 for 23 clus- 

ters, almost identical to the value obtained by Bahcall and 
Lubin (1994) from a larger compilation.) Clearly, either one 
or more of the assumptions of the above model is invalid or 
some of the measurements are in error. 

If we assume the "galaxies" to be distributed in the 
same way as the dark matter in our simulations, and we 
analyze them in the same way as the real data, we find a 
similar "/9-discrepancy" . When we fit the gas density profiles 
in the inner parts only, we obtain, ~ 0.6 (Figure 14), 
while a comparison of the velocity dispersion of the dark 



matter with the emission-weightedy^-Ta,y temperature yields 
/?T — 1 (column 2 of Table 3). This apparent discrepancy 
is readily explained as follows: (i) the model gas density 
profiles become steeper with radius so that if eqn (7) is fitted 
to the entire radial range (not just to the inner parts), the 
mean value of [ij for our six clusters is [ij = 0.82 ±0.06, not 
[ij = 0.6 (Table 2); (m) the dark matter density profile is 
shallower than King's approximation and has (/?/*) = 0.86± 
0.11 (Table 2) rather than [ij^ = 1. Hence Pf/Pf* ~ Pt 
expected. 

It is quite likely that a similar resolution of the "/?- 
discrepancy" applies to real clusters. As we showed in Sec- 
tion 4.1, the available data suggest that their gas density 
profiles differ from perfect /9-models in the same way as our 
simulations do. In addition, as Bahcall and Lubin (1994) 
have emphasized, the galaxy density profile in real clusters 
is better fit by ~ 0.8 than hy ^ 1. It is clear that 
since the gas and galaxy distributions in clusters are neither 
isothermal nor perfect fits to equation (7), an analysis in 
terms of these models must be treated with some caution. 

Our results are similar to those of Evrard (1990a, b), 
but our interpretation is slightly different. Evrard's clus- 
ters also exhibit a /9-discrepancy, but the values oi P ^ and 
Pj' he derives are significantly larger than ours. His value 
of Pj' ~ 1.2 exceeds ours because he used a mass- weighted 
temperature whereas we use an emission-weighted temper- 
ature. Since the models are not isothermal, the two weight- 
ings give different results. The emission-weighted tempera- 
ture is closer to the quantity measured for real clusters. In 
Table 3 we compare values of Pj' using emission-weighted 
(Tx) and mass- weighted (Tm) temperatures for the clusters 
as a whole (r20o) and for the inner regions alone (riooo)- 
(All columns in this table except the last one refer to the 
"equilibrium" configurations discussed in Section 4.2; riooo 
is the radius of a sphere of mean overdensity 1000.) Evrard 
ascribed the /9-discrepancy in his models to a combination 
of incomplete thermalization of the gas and to anisotropics 
in the velocity distribution of the dark matter. Incomplete 
thermalization is evident in our models also, but, as Table 3 
shows, it has a relatively weak effect on the estimate of Pj' 
when emission-weighted rather than mass-weighted temper- 
atures are used. Once these effects are taken into account, 
our results are in good agreement with Evrard's and with 
those of Thomas and Couchman (1992). 

The largest values of Pj' occur in clusters which have 
undergone recent mergers. Mergers stir up the gas in the 
outer parts and the resulting departure from hydrostatic 
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equilibrium can give rise to values of Pj' as large as ~ 1.4. 
This is illustrated in the last column of Table 3 which gives 
data for clusters identified at 2 = (within r20o)) rather 
than for the "equilibrium" configurations discussed above. 
At 2 = only CLl and CL6 are undisturbed and these 
have Pt{Tx) — 1- Provided the measured galaxy velocity 
dispersion is reliable, values of Pj' significantly in excess of 
one can be taken as an indicator of recent merger activity. 

None of our clusters has a Pj' value significantly below 
unity, yet there are several clusters in the samples of Edge 
and Stewart (1991a, b) and Lubin and Bahcall (1993), par- 
ticularly low-temperature clusters, with Pj' < 1. Unless the 
temperature of these clusters has been overestimated, such 
low values of Pj' indicate a bias in the galaxy velocity disper- 
sion similar to that seen in the simulations of Carlberg and 
Couchman (1989), Carlberg and Dubinski (1991), Couch- 
man and Carlberg (1992), Katz et al. (1992) and Evrard 
et al. (1994a). Such "velocity bias" should in principle be 
observable in the form of a steep number-density profile for 
the galaxies (i.e. in a large value of [ij^). 

5 "PREHEATED" GAS 

We have seen that a simple model in which the intr- 
acluster medium originates from an initially cold, nonra- 
diative gas leads to density and temperature profiles which 
resemble those observed in real clusters. Nevertheless, the 
X-ray luminosities of our clusters do not match observations 
very well. In particular, the predicted correlation between 
X-ray luminosity and temperature is too shallow {Lx oc 
rather than Lx'XjT^) and, for realistic cluster mass func- 
tions, the amplitude of the X-ray luminosity function in- 
creases with redshift, in conflict with observation. Since 
most of the X-ray luminosity originates in the cluster core, 
it seems worthwhile to investigate variants of the simple 
model which modify the state of the gas in the inner regions. 
One such variant, originally proposed by Kaiser (1991) and 
Evrard and Henry (1991), replaces the assumption that the 
gas is initially cold with the assumption that it was heated 
prior to cluster collapse by non-gravitational processes as- 
sociated with galaxy formation. Placed on a high initial 
adiabat, the gas in the central regions is less affected by 
shocks and evolves at nearly constant entropy. Eqns (1) 
and (3) give the scaling laws describing a cluster population 
in which all clusters are assumed to be isothermal and to 
have the same central specific entropy. 

To investigate this model in a quantitative fashion, we 
repeated the CLl, CL3 and CL5 simulations, this time set- 
ting the initial temperature of each gas particle so that its 
specific entropy was about the same as the final value at- 
tained in the core of CLl, our most massive cluster (see Fig- 
ure 4d). The resulting three simulations are labelled EPl, 
EP2 and EPS (see Table 1). The right-hand panel of Fig- 
ure 8 shows that the core regions of the preheated clusters 
do indeed evolve isentropically. As a result, the central gas 
density no longer exhibits the scaling behaviour of our ear- 
lier models. This is illustrated in Figure 15 where the thick 
solid line shows the ensemble fit to the density profile of the 
(initially cold) clusters of Figure 6b. The central overdensi- 
ties of the new models are lower and depend on the mass (or 
temperature) of the cluster. Cooler clusters have lower cen- 
tral overdensities than hotter clusters. Note, however, that 



changes due to the preheating are seen only in the central 
regions. In the outer cluster the density profiles are almost 
unchanged. 

Since the effects of preheating are greater in cooler clus- 
ters. X-ray luminosity increases more rapidly with mass and 
temperature in this model. This is shown in the two central 
panels of Figure 11 where the EP clusters are represented 
by solid triangles. The dotted lines indicate the prediction 
of eqn (3) for evolution at constant central entropy. The 
simulations agree extremely well with this prediction and 
give a much better fit to the observed Lx-T correlation 
(crosses and stars) than the initially cold clusters. Because 
of the numerical limitations discussed in Section 3.6, we can- 
not reliably determine the evolution of the X-ray luminosity 
function. However, the scaling relations obeyed by our clus- 
ters at 2 = imply that the Lx-T correlation should be 
independent of redshift and this, in turn, predicts a decline 
in the abundance of high luminosity clusters with increas- 
ing redshift. The rate of evolution depends on the initial 
fluctuation spectrum and for plausible choices it appears to 
be comparable to that seen in the limited datasets available 
(Evrard and Henry 1991, Castander et al. 1994). Finally, 
we note that the anticorrelation between the core radius (in 
units of r2oo) and mass seen in Figure 15 implies that the 
core radius of the gas distribution (in physical units) is es- 
sentially independent of mass. This agrees with the scaling, 
Tcore M~^l^' , expected if all clusters have the same cen- 
tral entropy (cf. eqns 1 and 3). 

6 CONCLUSIONS 

Our simulations of cluster formation make similar physi- 
cal assumptions and employ similar numerical techniques 
to those of Evrard (1990a, b). Nevertheless there are some 
significant differences; the present experiments have better 
resolution, follow clusters with a wider range of mass, and 
choose clusters according to the final nonlinear mass distri- 
bution rather than the initial linear overdensity field. Al- 
though our results agree well with Evrard's earlier work, 
we have chosen to address somewhat different issues. In 
particular, we have investigated whether, in the absence of 
cooling or heating, the distributions of mass and of gas can 
be treated as varying from cluster to cluster only through a 
scaling with size and age. We have also explored a model 
in which the gas is assumed to be heated before cluster col- 
lapse. Our main conclusions are the following. 

1) In a flat universe cluster formation by hierarchical 
clustering does indeed lead to objects whose properties scale 
at least approximately with size and age. Measuring radii 
in units of the radius enclosing a mean mass overdensity of 
200, densities in units of the cosmic mean, and temperatures 
in units of the corresponding virial value, we find that (to 
within a factor of 2 and over the limited radius range mod- 
elled reliably in our simulations) the density and tempera- 
ture profiles have a shape which depends at most weakly on 
time and on cluster mass. This result is expected, since the 
CDM power spectrum is virtually scale-free over the relevant 
range of length-scales and, in the absence of heating and 
cooling processes, there are no characteristic physical scales 
in the problem. Our success in verifying this scaling over a 
wide range of physical and numerical parameters suggests 
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that our techniques are adequate to resolve the structure of 
present-day clusters. 

2) Gas and dark matter evolve differently in the cores 
of clusters. Shocks cause a continual increase of gas entropy, 
while the maximum phase space density of the dark matter 
remains approximately constant. Numerical resolution is a 
significant limit on our ability to study dark matter cores but 
may not affect our estimates of the central thermodynamic 
properties of the gas. Although simulations with better res- 
olution are required to confirm this conjecture, we note that 
Pearce et al. (f994) come to very similar conclusions on the 
basis of their more schematic but higher resolution simula- 
tions of cluster merging. 

3) Our equilibrium clusters have X-ray properties which 
obey the standard scaling laws {e.g. Kaiser f986). In par- 
ticular, at 2 = our clusters satisfy Lx oc very well, 
and we find that clusters of a given mass are more luminous 
at higher redshift. Unfortunately, our simulations have in- 
sufficient resolution for a proper test of the scaling of lumi- 
nosity with redshift. Previous studies using Eulerian hydro- 
dynamic techniques (e.g. Kang et al. f994, Bryan et al. 
f 994) are likely to have predicted incorrect evolutionary be- 
haviour and an incorrect Lx - T relation because of poor 
resolution in the central regions of clusters. 

4) The dark matter and the gas have similar density 
profiles in our equilibrium clusters, but the gas is slightly 
more extended. The slope a = rflnp/rflnr of the gas den- 
sity profile increases gradually from a ~ — f at r ~ 0.05r2oo 
to a ~ —2.5 near r2oo- A standard /9-model yields [ij ~ 0.8 
if the whole range 0.03 < r/vjoo < f is used, but lower 
values are found if, as is common in observational studies, 
only the central regions are fitted. The gas temperature is 
almost constant in the central regions but tends to drop be- 
yond ~ 0.4r2oo- At T'joo it is typically half of the central 
value. These density and temperature profiles agree reason- 
ably well with those inferred for observed clusters without 
strong cooling flows. 

5) The parameter /Sy = /impa'^j^ /kT is close to unity 
for our equilibrium clusters provided emission- weighted tem- 
peratures are used. Larger values are often associated with 
ongoing mergers. None of our simulated clusters has Pj' 
much less than f, suggesting that "velocity bias" (i.e. a 
lower velocity dispersion for the galaxies than for the dark 
matter) may be substantial in those clusters observed to 
have P < 1. In our models, as in real data, Pj' > P^. As 
has been frequently noted (e.g. Sarazin f988, Evrard f990b, 
Bahcall and Lubin f 994) this infamous "/9-discrepancy" can 
be traced to the fact that the galaxy and dark matter dis- 
tributions are not well represented by King's approximation 
to the core of an isothermal sphere. 

6) Although formation models with a nonradiative and 
initially cold gas match most observed properties of clusters, 
their core structure results in X-ray luminosities which obey 
neither the observed scaling with temperature nor the ob- 
served evolutionary trend. As pointed out by Kaiser (f99f) 
and Evrard and Henry (f99f), models in which the gas ini- 
tially has high entropy are in much better agreement with 
observation. Such preheating only affects the core prop- 
erties of the gas; the outer density profiles are essentially 
unchanged. However, since observed cluster cores often 
have short cooling times, it is difficult to imagine that their 
specific entropy is constant and reflects processes occurring 



prior to cluster formation. It is unlikely that a full expla- 
nation of the X-ray properties of clusters will emerge before 
the role of cooling in cluster cores is properly understood. 

ACKNOWLEDGEMENTS 

We thank the referee, Peter Thomas, for a very careful read- 
ing of our paper. This work has been supported by the UK 
Science and Engineering Research Council. 

REFERENCES 

Arnaud, M., Rothenflug, R., Boulade, O., Vigroux, L., and 

Vagioni-Flam, E. 1992, Astr. Ap, 254 49. 
Bahcall, N., and Sarazin, C. 1977, Ap.J., 213, L99. 
Bahcall, N., and Lubin, L. 1994, Ap. J. Lett., in press. 
Bartlett, J.G., and Silk, J. 1993, Ap.L. Lett., 407, L45. 
Benz, W. 1990, in Numerical Modelling of Stellar Pulsations: 

Problems and Prospects, ed. J.R. Buchler, (Dordrecht: 

Kluwer) . 

Benz, W., Bowers, R.L., Cameron, A.G.W. and Press, W.H. 

1990, Ap.L., 348, 647. 
Bertschinger, E. 1985, Ap . L. Supp., 58, 39. 

Blanchard, A., Valls-Gabaud, D., and Mamon, G.A. 1992, Astr. 
Ap., 264 365. 

Bond, J.R., Cole, S., Efstathiou, G., and Kaiser, N. 1991, Ap.L., 
379, 440. 

Bower, R. 1991, M.N.R.A.S., 248, 332. 

Bower, R.G., BShringer, H., Briel, U.G., Elhs, R.S., Castander, 

F.J., and Couch, W.J. 1994, M.N.R.A.S., in press. 
Briel, U.G., Henry, J. P., and Bohringer, H. 1992, Astron. 

Astroph., 259, L31. 
Bryan, G.L., Cen, R., Norman, M.L., Ostriker, J. P., and Stone, 

J.M. 1994, Ap. L., in press. 
Carlberg, R.G., and Couchmann, H.M.P. 1989, Ap.L., 340, 47. 
Carlberg, R.G., Couchmann, H.M.P., and Thomas, P.A. 1990, 

Ap.L., 352, L29. 
Carlberg, R.G., and Dubinski, J., 1991, Ap.L., 369, 13. 
Castander, F.J., Ellis, R.S., Frenk, C.S., Dressier, A., and Gunn, 

J.E. 1994, Ap.L. Lett., in press. 
Cavaliere, A. and Fusco-Femiano, 1976, Astron. Astrophys., 49, 

137. 

Couchmann, H.M.P., and Carlberg, R.G. 1992, Ap.L., 389 453. 
David, L., Arnaud, K.A., Forman, W., and Jones, C. 1990, 
Ap.L., 356, 32. 

David, L., Slyz, A., Jones, C, Forman, W., Vrtilek, S.D., and 

Arnaud, K.A. 1993, Ap.L., 412, 479. 
Dubinski, J., and Carlberg, R.G. 1991, Ap.L., 378, 496. 
Edge, A.C., Stewart, 1991a, M.N.R.A.S., 252, 414. 
Edge, A.C., Stewart, 1991b, M.N.R.A.S., 252, 428. 
Edge, A.C., Stewart, G.C., Fabian, A.C., and Arnaud, K.A., 

1990, M.N.R.A.S., 245, 559. 
Efstathiou, G., Frenk, C, White, S.D.M. and Davis, M. 1988, 

M.N.R.A.S., 235, 715. 
Evrard, A.E. 1988, M.N.R.A.S., 235, 911. 

Evrard, A.E. 1990a, in Clusters of Galaxies, eds. M. Fitchett 
and W. Oegerle (Cambridge: Cambridge University Press), 
p. 287. 

Evrard, A.E. 1990b, Ap.L., 363, 349. 
Evrard, A.E. and Henry, J.P. 1991, Ap.L., 383, 95. 
Evrard, A.E., Summers, F., and Davis, M. 1994a, Ap .L., in 
press. 

Evrard, A.E., Mohr, J. J., Fabricant, D.G., and Geller, M.J. 

1994b, Ap.L. Lett, in press. 
Fabian, A.C., Canizares, C.R., and Nulsen, P.E.J. 1991, Astron. 

Ap. Rev., 2, 191. 
Fillmore, J.A. and Goldreich, P. 1984, Ap.L., 281, 1. 



ll-AUG-94 



12 Navarro^ Frenk and White 



Forman, W. and Jones, C. 1990, in Clusters of Galaxies, eds. 

M. Fitchett and W. Oegerle (Cambridge: Cambridge 

University Press) 
Frenk, C, White, S.D.M., Efstathiou, G., and Davis, M. 1990, 

Ap,J,, 351, 10. 

Gioia, LM., Henry, J. P., Maccacaro, T., Morris, S.L., Stocke, 
J.T., and Wolter, A. 1990, Ap,J, Lett,, 356, L35. 

Gunn, J.E. and Gott, J.R. 1972, Ap, J., 209, 1. 

Henry, J.P., and Arnaud K.A. 1991, Ap, J,, 372, 410. 

Henry, J. P., Gioia, I.M., Maccacaro, T., Morris, S.L., Stocke, 
J.T., and Wolter, A. 1992, Ap,J,, 386, 408. 

Hernquist, L. and Katz, N. 1989, Ap,J, Supp,, 70, 419. 

Hoffman, Y. 1988, Ap,J,, 329, 8. 

Hughes, J. P., Yamashita, K., Okumura, Y., Tsunemi, H., and 

Matsuoka, M. 1988a, Ap, J,, 327, 615. 
Hughes, J. P., Gorenstein, P., Fabricant, D. 1988b, Ap, J,, 329, 

82. 

Jones, C. and Forman, W. 1984, Ap,J,, 276, 38. 
Kaiser, N. 1986, M,N,R,A,S,, 222, 323. 
Kaiser, N. 1991, Ap,J,, 383, 104. 

Kang, H., Cen, R., Ostriker, J. P., and Ryu, D. 1994, Ap, J,, in 
press. 

Katz, N. and White, S.D.M. 1993, Ap,J,, 412, 455. 
Katz, N., Hernquist, L., and Weinberg, D. 1992, Ap, J, Lett, 
399, L109. 

Kauffmann, G. and White, S.D.M. 1993, M,N,R,A,S,, 261, 921. 
Lacey, C. & Cole, S. 1993, M,N,R,A,S,, 262, 627. 
Lubin, L.M, and Bahcaff, N.A. 1993, Ap,J, Lett,, 415, L17. 
Monaghan, J.J. 1992, Ann, Rev, Astr, Ap,, 30, 543. 
Mulchaey, J.S., Davis, D.S., Mushotzky, R.F., and Burstein, D. 

1993, Ap,J,, 404, L9. 
Mushotzky, R. 1984, Rhys, Scripta, T7, 157. 
Navarro, J.F. and White, S.D.M. 1993, M,N,R,A,S,, 265, 271. 
Ponman, T.J., and Bertram, D. 1993, Nature, 363, 51. 
Quinn, P.J., Salmon, J.K, and Zurek, W.H. 1986, Nature, 322, 

329. 

Pearce, F., Thomas, P.A., and Couchman, H.M.P. 1994, 

M,N,R,A,S,, 268, 953. 
Press, W.H., and Schechter, P. 1974, Ap,J,, 187, 425. 
Richstone, D., Loeb, A., and Turner, E.L. 1992, Ap,J,, 393 477. 
Sarazin, C, 1988, in X-ray emissions from clusters of galaxies, 

(Cambridge:Cambridge University Press). 
Spitzer, L. 1978, Diffuse Matter in Space, (New York: Wiley). 
Thomas, P.A., and Couchman, H.M.P. 1992, M,N,R,A,S,, 257, 

11. 

Watt, M.P., Ponman, T.J., Bertram, D., Eyles, C.J., Skinner, 
G.K., and Willmore, A.P. 1992, Mon, Not, R, ast, Soc, 
258, 738. 

White, S.D.M. 1982, in Morphology and Dynamics of Galaxies, 
Saas-Fee Lectures, ed. L. Martinet and M. Mayor, 
(Saas-Fee:Geneva Observatory). 

White, S.D.M., Navarro, J.F., Evrard, A.E., and Frenk, C.S. 
1993, Nature, 366, 429. 

This paper has been produced using the Blackwell 
Scientific Publications macros. 



Figure 1: Evolution of the distribution of gas (left) and dark 
matter (right) for cluster CLl. This is the most massive cluster 
in our ensemble of simulations. The epochs shown correspond to 
redshifts 2, 1, 0.4, and 0. Units are in physical kpc for 
i^Q = 50 km s~^ Mpc~^. Note the resemblance between the 
distributions of gas and dark matter at all times. 



Figure 2: The gas distribution at the final time i^z = 0) in 
clusters from the CL series. The dark matter distribution (not 
shown) is similar to that of the gas in all cases. Units are in 
physical kpc for Hq = 50 km s~^ Mpc~^. Note that 4 out of 
the Q CL clusters are undergoing a major merger event at 2: = 0. 



Figure 3: The evolution of the dark mass inside a sphere of 
current mean overdensity 200 centred on the most massive 
clump. Masses are in units of the overdensity 200 mass at 2: = 0. 
Most of the cluster mass is accreted during merger events, 
reflected in sudden jumps in mass. Clusters form late; three of 
them assemble half of their mass only after z = 0.2. 



Figure 4: Evolution of the hydrodynamic properties of the 
intracluster gas in the six CL clusters. The "density" and 
"temperature" panels show the time evolution of the core (solid 
line), main body (dotted line), and outer boundary (dashed line) 
of the cluster. These properties were obtained by averaging over 
10% of the particles which are nearest the core, half-mass radius 
and outer boundary respectively. The gas density and entropy 
are in code units (§2.3). For clarity, the temperature curves have 
been shifted downwards successively by one order of magnitude. 
The X-ray luminosity includes the emission due to all particles 
which at 2 = make up each cluster. The "entropy" panel 
shows the evolution of the core and boundary of each system 
only. The initial entropy of "pre-heated" clusters i^EP series) is 
shown by a horizontal dotted line. Units are given in §2.3. 



Figure 5: Gas and dark matter density profiles for equilibrium 
configurations oi CL clusters. The profiles were constructed by 
averaging the density in spherical shells containing 50 particles 
each. The units of density are given in §2.3. The different curves 
show different clusters in order of decreasing mass. 
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Figure 6: Scaled density profiles of CL clusters. The solid line 
segment in the upper panel has the slope, p oc of an 

isothermal sphere. Gravitational softenings are shown as small 
vertical lines with the same line type as the corresponding 
profile. The lower panel also shows the best ensemble fit using 
the isothermal /5-model (thick solid line) and Bertschinger's 
(1985) self-similar solution (thick dotted line). 



Figure 7: As Figure 6, but for clusters identified at 2: = 1. The 
thick curves are the same as those shown in Figure 6. Note that 
at this earlier time, the profiles deviate from the z = ensemble 
fits at increasingly large radii because the number of particles in 
each cluster is smaller and the gravitational softening lengths 
(shown as vertical lines) become an increasingly large fraction of 
7^200- Outside the core, the profiles are very similar to those at 
z = 0. 



Figure 8: Evolution of core particles (innermost 10% of 
the mass at 2: = 0) in the density-temperature (gas) and 
density- velocity dispersion (DM) plane. Dotted lines show loci 
of constant gas entropy, T / p^/^, and constant dark matter 
phase-space density, p/<7^. Left and right panels show the 
evolution between z = 1 (open symbols) and z = (filled 
symbols) of CL and EP clusters, respectively. Units are 
arbitrary and relative values of the "entropy" increase upwards. 



Figure 9: Scaled temperature profiles of the CL clusters. The 
thick dotted line shows Bertshinger's (1985) self-similar solution. 
The thick solid segment indicates T oc Vertical segments 

show the gravitational softening. All clusters are nearly 
isothermal out to r ~ 0.4^200 and typical temperatures at the 
virial radius are about one-half of the values near the center. 
The units are described in the text. 



Figure 10: The ratio of binding mass (estimated from eqn 9) 
to total mass as a function of radius. Incomplete thermalization 
leads to a systematic underestimate of the true total mass. 



Figure 11: Correlations between the bolometric X-ray 
luminosity, total mass (at mean density contrast 200), 1-D dark 
matter velocity dispersion, and X-ray emission-weighted 
temperature, for CL clusters (filled circles) and pre-heated 
clusters (filled triangles). Observational data compiled from 
sources quoted in the text are plotted as small crosses and 
starred polygons. Solid lines show the scaling relations (eqns 1 
and 2) appropriate to an initially cold gas; dotted lines show the 
corresponding relations (eqns 1 and 3) for the case when all 
clusters have the same initial minimum entropy. Open circles in 
the luminosity-mass panel represent the progenitors of CL 
clusters at 2 = 1 and the dashed line shows the evolution 
expected from the scaling relations. The weaker evolution seen 
in the simulations reflects the fact that clusters identified at 
high redshift are more poorly resolved. 



Figure 12: Scaled density and temperature profiles for 
simulations with the same initial conditions as CLl but with 3 
times the total number of particles (70000 particles in total; 
dotted line) and 1/2 the total number of particles (11000 
particles in total; dashed line). The solid line shows the profiles 
for the CLl cluster. No significant differences are noticeable. 



Figure 13: The dependence of the X-ray luminosity of cluster 
CLl at 2 = on the number of particles used in the simulation. 
Luminosities are in erg/s, and Ngas is the number of gas 
particles inside the virial radius of the cluster. In order of 
decreasing Ngas, the points correspond to runs CLl X 3, CLl, 
CLl/2, (see Table 1). The other two runs correspond to 
decreasing the number of particles of the CLl run by factors of 
3 and 6, respectively. 



Figure 14: The dependence of the parameter j3j (eqn 7) on 
the outermost radius, rmax, used in the fit. Open circles 
represent the simulated clusters. Filled circles correspond to the 
eight clusters from Jones and Forman (1984) for which rmax 
can be estimated (see text and their Figures 1 and 2) and to the 
Coma cluster from the results of Brief et al. (1992). In both 
models and real data, lower values ofj3j are obtained for lower 
values of rmax- 
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Figure 15: As Figure 6, but for clusters with pre-heated gas 
(EP in Table 1). The three clusters have the same initial 
central entropy. More massive clusters have higher central 
densities and smaller core radii (in units of their virial radii) 
than smaller clusters. 
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